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applicability of most numerical techniques including the 
Finite Element Method, when very accurate results are required. 
This is particularly true when dealing with Boundary Value 
Problems. In order to surpass this difficulty a new method 
to solve these problems has been devised which does not re- 
quire matrix storage while still providing the possibility of 
accuracy improvement. 
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lL. SON ERODUCTION 


For many engineering problems it is not always possible 
to find an exact solution. An exact solution is an analyti- 
cal mathematical expression that gives the value of the unknown 
at any point within a previously specified range. For problems 
involving complex material properties and boundary conditions, 
different numerical techniques have been developed to approxi- 
mate the exact solution to a more or less acceptable degree 
of accuracy. One of these techniques is the Finite Element 
Method, where a given problem is to be discretized, that is to 
say, approximate solutions are to be found at discrete points 
in the body. The way one accomplishes this, is by subdividing 
the whole body into finite elements. The solution is then 
formulated for each small unit and combined to obtain the 
solution of the whole system. Obviously the greater the number 
of elements, the better the accuracy. 

However, despite the fact that high-speed digital computers 
have enabled engineers to successfully apply these numerical 
techniques, we still face the problem of large matrix storage. 
Complex problems requiring very accurate approximate solutions 
lead to large matrices and consequently larger computer memory 
will be required. This fact constitutes a limitation on the 


applicability of these numerical techniques. 


The method we are dealing with in this work intends to 
EN cEthossdisrtficultyset гейпсіпа а В.У.Р. to ап І.У.Р. Unlike 
the "shooting" method, where basically the same idea is used 
together with an iterative scheme to achieve a solution, the 
method developed here achieves a solution without iteration. 

Although restricted to a particular type of linear differ- 
ential equations and only to one-dimensional problems, as it 
will be shown, this method can be applied to many problems 
where no exact integrable forcing functions occur, or, even 
more, the function exists only as a set of values at specific 
points. 

Basically what the method does is approximate the exact 
solution by a set of linear functions, each of them applying 
in a very small interval. However, these functions are not 
independent from each other. We use the previous one to find 
the next, until we eventually reach the other end of the range 
in consideration. Depending on the order of the differential 
equation, we will find as many sets of linear functions as 
required by its order, each set applying to a specific inte- 
gration, and as before, every set depends on the previous one. 
The way in which we will do this will allow us to correlate 
all possible boundary conditions in an explicit manner which 
transforms boundary value problems into initial value problems 
and by doing that, get to the result. 

The method is indeed an approximation and as such, it is 


subject to error. We can say, however, that by decreasing the 


step size we get better accuracy, provided we carry enough 
Significant digits to take care of round-off ево АЕН 
respect to the later a good programming technigue is required: 
Finally we shall apply the method to several practical applica- 
tions, specifically to beam problems where a fourth order 
differential equation occurs, and a variety of boundary 


conditions can be given. 


II. BACKGROUND 


A. GENERAL 
In this section we are going to introduce the basic ideas 


of the method we will be dealing with. Let 


Ү'(х) = f(x) (21) 


(еге 15 understood that the function f(x) may or may not 
Se Exactly integrable, but we can integrate it numerically. 


From the above relationship we have: 


ИСКИ E Y (252) 


Here the constant of integration Y is assumed to be zero 
for the moment. In the most general case, it is clear that 
there is no way to know what the function Y(x) is, since we 
may be dealing with non-integrable functions. However, we 
can approximate the function Y(x) by a linear function y(x), 
provided the interval in which this approximation applies is 


sufficiently small. So we take 


Yoo. — a yak) = ес (298) 


Now, by following the general procedure of integration, and 


from (2.2), we have: 


b 
юх + с| = fJ f(x)dx (2.4) 


where (a,b) denote the limits of the interval under considera- 


tion. After simplification Irt Is round that: 


Кт хах 
m = 一 一 一 一 一 一 2-5) 


b -a 
that is to say, the slope of the approximating line given by 
(2.3) can be explicitly determined. Actually, as we will see 
later, it can be said that the "general" approximate solution 
to (2.1) has been found within the interval from a to b. In 
order to determine the "particular" solution, in other words 
to find the intercept c of (2.3), an auxiliary condition must 


be imposed. Let 
Y(a) = Yo 
then: 


У (а) = У = ma с 


or 


(210 


|^ 
X 
| 入 
oO 


y (x) m(x-a) * Yg a 


l0 


where m is given by (2.4). Now we have fully determined the 
approximate solution inside the given interval. Note that 
(2.6) will give results at points a and b, as accurate as the 
numerical integration performed on (2.5). 

Fig. 2.1 shows the whole process so far. The upper curve 
represents f(x), which is given. The lower curve is Y(x), 
the integral of f(x); however within the interval a to b, Y(x) 
1S approximated by a line. Point Yo is the given auxiliary 
condition and point Y, Can be determined exactly from (2.6), 
by letting x = b. 

The next step now becomes evident, since point Үр сап be 
determined. We are now in the same position as before, so 
all we need to do is repeat the process over again. However, 
this time with a new auxiliary condition; the last one we 
have just found, and over the next interval. We keep going 
this way until we reach the other extreme of the range where 
the differential equation applies. See Fig. 2.2. 

In summary, we can say that every step we take we are 
solving a "new" differential equationby approximating lines 
which gives us true values at the points of intersection. 

The fact that we have a "new" differential equation at every 
single step allows us to deal with discontinuous functions and/ 


or functions existing only as a set of values. 


EC THE FIRST INTEGRATION 


AE DIS Point let US Introduce a new parameter: 


h = b-a (2.7) 


ы) 


o< 


C ME 


FIGURE 2.1 
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f (x) 


o< 





FIGURE 2.2 
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where h is assumed to be small and constant throughout the 
analysis. We are now in a position to determine a general 
relationship which will allow us to find the value of Y(x) at 


discrete points, namely at the intersection of the strada 


lines. From (2.6) we have: 
У (ВА - пра Y. = mh + ze. 
In general, 
xm = m.h + Yi 0 = Dea (2990 
where 
xX. 
1 
f f(x)dx 
X. 
2 у 
m; = Б (259) 
Би: 
E S n 4 "аа 
- = m. oP + Yi. 3 
У = mh + Yo 


Replacing these last relationships in (2.8) and after factor- 


ing h, we have: 
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Or 
) 
Ү - h т. + Y 210) 
n gary О 
Emwlarly, 
h 
POS: 
m = 9 
1 h 
2h 
1 cres 
|. di 
i ne x h 
Adding terms gives:: 
nh 
в a Е (9 а= 
р D s (2L) 
1=1 


Dicam replace in (2.10) to obtain: 


EE ВЕ (212) 


BMMMCE X = nh. Equations (2.10) and (2.12) constitute two very 
important results. Eq. (2.12) is indeed not an unexpected one, 
and several considerations can be drawn from it. First of all 


the value of "m does not depend on the step size h, at least 
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for this first integration, so its accuracy is determined only 
by the exactness of the integration routine used. Second- 
values of Y(x) can be directly determine at any arbitrary Tonki: 
n, that is to say we need not find its previous values as is 
the case for many numerical methods. Later this fact will 
become more useful. Finally, we see that, by transposing terms 
we can find either the left or right B.C., and Ећа5 арравеш 
unimportant fact is indeed a key step, since as we will see 
later, we will be able to correlate B.C.'s of higher order 
differential equations and solve for them without having to 
generate all previous values of the unknown function. In other 
words, we could transform a Boundary Value Problem into an 
Initial Value Problem and vice versa, depending on what is 


known, and what we are looking for. 
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III. HIGHER ORDER DIFFERENTIAL EQUATIONS 


A. GENERAL 

The previous chapter was dealing basically with the main 
ideas of the method, and we have solved a first order differ- 
ential equation. We are now going to extend the method to 
higher order equations. Essentially what we will do is, inte- 
grate several times the right hand side of the given equation; 
as required for its degree, using the same approach as before. 
However, we should keep in mind that it is only for the first 
integration that we will use the given function f(x). After 
this integration, this function is no longer available because 
what we have is a set of straight lines, each applying to a 
Specific interval. We need to find now a method that will 
allow us to perform a second integration of the function f(x), 
from the given set of lines. 

In order to do this, let us consider the curves of Fig. 
3.1. Suppose we have been given the following second order 


differential equation: 


Y"(x) = f, (x) (3:1) 
represented in Fig. 3.1 by the upper curve. Let us assume we 
have determined Y'(x) = f} (x) after a first integration by the 


method just introduced, Let 
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f (x) 
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X 
УЬ 
X, 
b 
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gy 
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uc) m £x) 


be the exact result, which is the lower curve in the same 


figure. Now let 


y'(x) = тух + Cy 
and 

ух) = Mx TE 
be the approximating lines to £1 (x) and f(x), respectively, 
that applies inside the interval from a to b only. In this 


figure, we can easily see that the area enclosed by the upper 
straight line is almost equal to the area enclosed by the 
lunction Е. (х). Obviously as the interval becomes smaller and 
smaller, both areas tend to be equal; in the limit they are 
indeed equal. It is evident now that we can approximate the 
exact integration of f (x), by integrating the line 


y'(x) = mx + cC provided the interval is small enough so the 


gu 
error is negligible. Now, as before, we need to determine m 
and c. Recall that m and c, have already been found. So by 


virtue of (2.5): 
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Бик рулы 2.1) 


b 
| (тух + с))ах 
m = Е 
ог: 
Rl 
m = =z (а == Ы) + C1 (3.2 
So we know now the slope of the next line; the line belong- 
ing to the solution function. In order to find the interco 


c, another auxiliary condition must be supplied. Let: 


Y(a) = Уа 

then, replacing this condition in the line y(x) = mx + с, ме 
Cem. 

V - ал = па + с 

а 

It can be shown that: 

у(х) = mixes Y (Зе 
where m is given by (3.2). Eq. (3.3) gives the approximate 


solution of the original second order differential equation, 


inside the interval a to b. Note that both auxiliary conditions 
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have been given at the left hand side. However, it really 
need not be like this, as will be shown later. 

At this point some considerations are in order. With 
reference to Fig. 3.1 where the upper curve is the given 
munction Y"(x) - f(x), the middle curve is Y'(x) = f (x), 
and the lower one is the solution Y(x) = f(x). 

As explained previously we have approximated the solution 
by straight lines. Consider now the slope of the lower 
ELErarsght line myghich is given by (3.2). This equation can be 


written as: 


= sl (m,a + с.) + (mb + са)! 
П 
mM = zly (a) 本 y, (5)] 
El 
I JAM r Yp] СОА) 


Since the ordinate at any point x in the Y'(x) curve is 
indeed the slope of the Y(x) curve at the same point, then 
by virtue of (3.4), we can say that the slope m of the lower 
straight line is the average of the slopes corresponding to 
its end points, namely point Y. and Y, тит 3.1. Recall 
that these two slopes are given exactly by the ordinates of 
the middle curve Y'(x). Furthermore let us evaluate 
y'(x) = тух + с), the equation of the upper straight line, at 


the midpoint of the interval, that is we let: 


21)! 








X = = (а * b) 
in the line y^ (x) and we get: 
十 
vě D T 4 Е 2 m 


As seen, it reproduces m, the slope of the lower straight 
line. 

Now we can have a deeper insight of the whole process, 
and why it works. See that what we are really going to do is 
to use the slope at the midpoint (the average of the slopes at 
the extremes) to approximate the true slope of the straight 
line y(x), and later determine the corresponding intercepts. 
In fact, had we known exactly the function Y'(x) = E we 
would have been able to determine the true slope of v(x) by 
direct integration as before. As we can see now, we are intro- 
ducing a source of error, and from now on we will not have as 


accurate results as in the first integration. 


B. DEVELOPMENTS OF TEE Mob 

Before we go into a general development of this method a 
slightly different notation must be introduced, and we will 
be using it throughout the rest of this research. Tne type 
of differential equations we will be looking at are of the 


form: 


ү М ху = f(x) (575 
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Such equations arise in the study of deflections of beams, 
and we shall apply the method to practical beam problems. Let 
us introduce the notation to be followed. We will write the 


approximate solution line of (3.5) as: 


у(х) = mx + C 


in a clean way without subscripts. The approximation line 


to the first derivative will then be: 


у'(х) = тух + С, 
and for the second derivative: 
y"(x) = m4X + с. 


and so on. 


The boundary conditions will be represented as: 


У (Хо) = Y Y(x,) и Y, 
Y' (x) E Y5 77 (554 = ps 
Y" (x) = YA ей 22-7 = ЩЕ 


etc., where the subscript o refers to the left hand side, and 


the subscript n to the right, of the range under consideration. 
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There is another subscript that we must keep track of. 
Recall that the whole solution function is going to be repre- 
sented by a set of approximating lines, see Fig. 2.2: Then if 
we start at the left end, we shall use the subscript n as a 
second indice to identify each line as we move rightwards. Be 
aware that we will take the subscript i immediately to the 
right of each line to name it locally. 

Some words must be said about the auxiliary conditions. 
Recall that in a specific problem the number of these condi- 
tions are equal to the order of the differential equation to 
be solved. If all these conditions are given at the same 
point, then we are dealing with an Initial Value Problem. 
However, if the auxiliary conditions are given at both ends 
of the range of interest, then we have a Boundary Value Prob- 
lem. When we are able to have an exact solution, we usually 
are led to a set of simultaneous equations involving the unknown 
B.C.'s; those can be solved algebraically and lead to the 
whole solution. However for equations with no exact solution, 
like the ones we are interested in, there is no way to deter- 
mine the conditions at the other end of the interval, and that 
is precisely the goal of our approach. We need to correlate 
in an explicit way the known conditions and the unknown ones 
by means of equations of the form of (2.10), so we can be able 
to find them. 

Note that in this equation the term h | m. is independem 


ісі і 
Of B.C.'s and it can be readily found even if we do not know 
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them. By similarity to an exact problem we can look at this 
term as being the general solution of the differential equation. 
Note that this term actually links the left and the right 
boundary conditions in a very explicit way. This is the 

heart of what the present method is all about. With this 

idea in mind we are now going to find several relationships 

for successive integrations. 

Let us start by noticing that (2.10) is a general relation 
that applies not only to the first integration but also to 
successive ones. The only term that is particular for each 
integration is the summation of the slopes. So in order to 
perform the next integration we shall determine this term 
specifically and substitute it into (2.10). We are now dealing 
with a fourth order differential equation since this type 
occurs in beam problems. The application to other orders will 


be self-evident. Let 
ОО = f(x) (3:6) 


After a first integration we have: 


n 
"n 一 ит з 7 
ще h Да rods (S) 
i=l 
Ducere the subscript 3 in the summation indicates that these 
slopes belong to ¥" = £5(x) and are given by (2.11). We need 
to find now 
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n 
Y mz Do ШЕ. (3-0 
і=1 


Following our convention, m is given by: 





2,i 
ш 1 
Dod. wm ЖО с. 
but 
= "1 Е 
be NEN M И 
then 
uS 
7 "1 m 
m = 2 m 
ог 
pos 
са 7 , "i 
m5 i — 5 ћ + vil (3.9) 
Now, since 
п- 1 
"1 = "1 
Үп = ШЕЕ 
des 


then 
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De п-і 


)a п, + mu (З 0)‏ 2 ي 


і=1 


This last relationship in the way it has been derived 
becomes a general one, and expresses the slope m n at any 
7 
point i as a function of the previous slopes. Similar rela- 


tions can be derived for successive integrations just by 


shifting the first indice. Now we need to find the summation 
of these slopes. It becomes: 
pem n-l ji 
уп = 5 ЈЕ m. *h ) lc т. E Gel 
ізі 1=1 a 


New, replacing (3.11) in (3.8): 


h n n-l 1i 
E в [5 3m th 2 | m. jes ] * Ys 
gb - 
1=1 1-1 
Or 
2.1 n n=l ji 
" = Ll ча 3 3.12 
h [5 ) 5 m. + L ) з т] pu LEY. ( ) 
Pel 
1=1 ізі 
Let: 
1 n n-l 1 
P та, |5 
1=1 i=l j=l 
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Note that, as stated above, this last term 52 15 also 
independent of boundary conditions and can be readily deter- 
mined. Equation (3.12) gives the particular solution toma 
second order differential equation as well as, in this 
specific case, it represents the second derivative of Y(x). 

The whole process, so far, can be summarized as follows: 
Equation (2.10) gives the solution at any point n of any 
integration, first second, etc. In order to use this equation 
we need to find the summation term, and this is given by a 
relation of the form of (3.11), which is also a general rela- 
tionship and expresses the summation of slopes in terms of 
the previous one. Recall that the only slope summation we know 
is the one involving the known function f(x), and is given by 
(2.11). That is why all further summations must be expressed 
as a function of this one. This last summation is then re- 
placed in (3.8), and simplified if possible. The next step is 
to identify the terms that are independent of boundary condi- 
tions and isolate them as in (3.12) and to evalute them 
separately. 

We can continue in the same way until we get a solution 
formula for Y(x). It can be shown that the complete solution 
for a fourth order differential equation of the form v(x) = £(x) 
is given by the following relationships: 


p = 5. + Ше (3.133 


9 


" = из " 
Үр = hs. + ППУ. + Үс (3.14) 


2Е 


where: 


and 
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+ Nh 
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1 ј=1 к=1 


те = к= E 
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(3486) 


ИК 
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(52-221), 


(3322) 
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IV. ERROR ANALYSIS AND CONVERGENCE 


Refer to Fig. 4.1 where we have plotted the functions 
Ү' = f- (x) and Y = f(x) together with their respective 
approximating straight lines inside the interval x = a to 
x = b. We assume that the upper straight line belongs to a 
first integration; in other words we have not introduced any 
error so far. A second integration, however, will be performed 
using this line y4 (x) since we do not know the function f- (x). 
It is worth considering now what happens when we integrate 


the line instead of the function itself. Let m be the slope 


of the lower straight line. This is given by: 


b 
J £, (x) dx 


= a 
b-a 


Let m* be the approximating slope given by: 


ES M ee ыс (4.1) 


fo (myx + c,)dx > | Е (ах (4.2) 
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since the area enclosed by the curve is smaller than the area 
enclosed by the line (we assume that there are no inflection 
points and/or discontinuities inside the interval). It is 


now clear that: 


m > m 
and the error introduced is: 


b b 
ӨШ Buc е с соја (1) 
a-b a a 


which is egual to the shaded area in Fig. 4.l. 

As indicated in Chapter III, m* comes out to be the average 
of the slopes at the end points of the lower straight line, 
and it is this slope we work with. So, the actual approximat- 
ing line y*(x) = m*x + c*, shown in Fig. 4.2 is steeper than 
y(x) since m* > m. The distance ED becomes the error introduced 
and we can not evaluate it because we do not know f(x). Another 
integration using y*(x) will make things even worse. Note 
that the error to be committed this time will be larger and 
is given by the shaded area in Fig. 4.2. 

So far it appears that the method is divergent in nature 
due to the fact that the error will grow bigger and bigger 
unless we find some means to correct it. At this point, the 


only possible way to do this is by reducing the "shaded" 
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areas, which means that we have to reduce the step size. 
Fortunately, we can say, at least in theory, that in the limit 
there will be no error; in practice however, a very small 
Step size will increase the round-off errors. The examples 
provided show that the approximation is indeed acceptable. 
There is a special case, however, where an error correc- 
tion can be made. This applies only to boundary value problems. 
Recall that a solution of any differential equation by this 


method is given by the general relationship: 


In B.V.P., we know the conditions at both ends of the interval. 


So this last equation can be written as: 
Е О n 
) m, = ——— (4.4) 


where the right hand side is known. Now since we are actually 
dealing with approximations to the true slopes, what we really 


fave as a slimmation is the next term, Call it M: 
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Our approximate solution is then: 


n 
* = 
үз DM ee (4.5 
i-1 


and the error at point n will be given by: 
e = ух - у (4.6) 
n 
Now let us see how the correction could be carried out 
individually at each interval. To do this we must find the 


exact slope m. It is given by: 


m = мх = e 


where e is the shaded area of Fig. 4.l and is given by (4.3). 


The slope summation is therefore: 


or 
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Хан — jhe (4:7) 


J J 
When j = n, we have: 
= х= 
ха Y nhe 
Or 
У = nhe 
n n 


but from (4.6) we get: 


= „A 
ез - ЗА (4.8) 


and the general relationship given by (2.10) transforms into: 


j | 
Y. = по) му - Је (4.9) 
j iE 1 o mm 
qd 
which is the corrected solution at any point j. In summary, 
what we need to do to correct the solution is: a) determine 


Y by performing the integration without any correction by 
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using equation (4.5); b) compare Y withitne given B C: ҮМ 
and obtain e; c) Perform again the integration this timed EM 
the relationship given by (4.9). 

From (4.9) we can find a whole set of corrected equations 
similar to those given by (3.13) to (3.20). However, Dy dio EEI 
this we will enormously complicate those relationships and 
the computational effort, together with round-oft errors ne 
not give any advantage at all, especially since we expect the 
local error e to be very small. 

There is, however, another stronger reason not to do that. 
The fact is that we have assumed that the error e is a con- 
stant and this is not always the case. If we restrict our- 
selves to integrate linear and/or constant functions, then the 
corrected method could be justified since for these functions 
the error given by (4.3) is a constant, for the shaded area 
of Fig. 4.1 will always be the same. In more general cases, 
the error e will be unpredictable and no correction can be 
made. 

This example serves only to illustrate how the error 
correction could be carried out. But it would not be practical 
and a reduction in step size appears to be the best correction 
as we will see in the examples. 

As stated above, the method appears to be divergent in 
nature, so our approximate solution will look like Fig. 4.3. 
In this figure, note that for any integration after the first 


one, the set of approximating lines diverge from the exact 
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solution, the further we move to the right, the greater the 
error introduced. In dealing with I.V.P.'s, the points Ye 
and Y, are known and these are our starting conditions. The 
solution will appear then as shown in this figure. Since we 
do not know the exact end points at the right hand side, 
namely points Y' and Y, there remains an uncertainty in the 
accuracy of the solution. 

For B.V.P.'s, however, the other extreme points are known. 
Referring to Fig. 4.3, suppose we are given as boundary condi- 
tions, points x and Y in the lower curve of this figure. 

The method we are dealing with requires that we know the 
initial points to be able to start the integrating algorithm. 
If we know the exact starting points, namely Mé and Yo’ then 

by using these initial conditions, our solution will appear 

as shown in Fig. 4.3. That is to say we will not end up at 
point Y which is exact, but at point Y*. То be able to reach 
the exact point Y, which we assume is known, we need to give 
the algorithm a "wrong" starting point E It is this appro» 
mate starting point which we find using the relationships given 
by equations (3.13) to (3.22) that allow us to match the exact 
solution in the whole range as shown in Fig. 4.4 in the lower 
curve. 

Note that if we were using corrected relationships we 
would be able to supply the algorithm with "exact" starting 
values, and still match the exact end points. Unfortunately, 


it is now the upper set of solution lines that absorb the 
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inherent error of the method since we have given them approxi- 
mate initial values. The solution now is going to look like 
Fig. 4.4. The lower solution becomes "exact" but the upper 
actually shifts apart from the exact. We can explain this by 
saying that the known B.C.'s are fixed by the algorithm while 
the free (unknown) ones will take the error, and this is 
exactly what happens in the examples in the next section. In 
summary, all we do is to shift the error from one place to 
another. We can not get rid of it unless we use corrected 
relationships or an infinite number of intervals to minimize 


the local error. 
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V. APPLICATIONS AND RESULTS 


A. GENERAL 

In this chapter we are going to study a particular appli- 
cation of the integrating algorithm introduced in the previous 
sections, namely, the deflection of beams, where a fourth 
order linear differential equation occurs, and several combina- 
tions of B.C.'s can be considered. 

From Mechanics of Solids, the relationships governing the 


deflection of beams are given by: 


EIY "(x) = а(х) 人 (5.1) 
EIY" (x) = V(x) (shear force) @ 
EIY"(x) = М(х) (moment) 9) 
Y'(x) = slope 
y(x) = deflection 


The method we are dealing with is basically a method of 
successive integration. Because of that we can start to inte- 
grate from any of the above relationships provided we know 


explicitly the right hand side of the equation. In some cases, 
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a given problem could be solved by integrating the moment 
second order differential equation. By doing this, we would 
be able to achieve more accuracy since only two integrations 
need to be performed and the accumulated error will be rela- 
tively small. 

However, when the loading of the beam is a complicated 
distribution and the expression for the bending moment is 
difficult to obtain, then the fourth order differential equa- 
tion will be the one to use. One advantage in using this 
expression is that the integrating algorithm will provide 
results for the shear force, bending moment, slope and 
deflection simultaneously. When using the second order equa- 
tion we only get slope and deflection. 

Since eq. (5.1) is a more general relationship, we are 
going to use this expression in the examples. In fact, most 
problems can be expressed in that way. Recall that a concen- 
trated force can be treated as a distributed force acting in 
a very small interval, and a concentrated moment reduces to a 
couple of concentrated loads acting in opposite directions a 
small distance apart. Example 2 illustrates this procedure. 

The examples provided in this chapter account for all 
types of B.C.'s. Fig. 5.1 shows the 4 cases to be treated, 
and we shall study them one by one in the next subsections. 
In all cases the flexural rigidity EI and the length L is 
set equal to l, for simplicity. Discontinuous types of loading 
are emphasized, and of course, many combinations of loadings 


can be dealt with by simply using the principle of superposition. 
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The computer program used to solve these problems is shown 
at the end of this chapter, and it is a straightforward cock- 
of equations (3.13) to (3.22). This programming technique ни 
particular, tries to minimize round-off errors. For now, it 
suffices to say that all problems are solved using SINGLE 
PRECISION, one hundred intervals (problem 2 uses 200), and we 
get accuracy to the fourth and even to the fifth decimal places 
in some examples. Finally, only the fourth integration, the 
deflection, is compared with the exact solution. Recall that 
the error is cumulative, and it is here that we expect the 
bigger error. We may conclude then that preceding integrations 
are more exact. However, we should keep in mind that an 
approximated starting point is needed in some cases, and the 
results will be shifted by some amount from the exact. We 


will see this as we proceed into the next section. 


B. EXAMPLES 
1. Example 1; @¥ameed=lrecemecam 
Refer to Fig. 5.la. This is a statically determina е 
cantilever beam loaded over half its length by a uniformly 
distributed load of 1 unit of weight per unit of length as 
shown. We shall determine the shearing force, bending moment, 
slope and deflection at discrete points of the beam. 
2 поста тош 
This is a B.V.P. but for purpose of illustrating how 
to use the method when dealing with I.V.P.'s, we will treat 


it as such. All the required initial conditions can be drawn 
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from the given configuration. Here the function f(x) to be 
integrated is discontinuous and given by a set of two differ- 


ent functions defined by: 
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We need to determine the initial values to start the 
algorithm. From statics, it can easily be shown that the 


mittital conditions are as follows: 


Shear: Y = -0.5 ша = 0.0 
Moment: Y = 0.375 Үр = 0.0 
Slope: УС = 0.0 ү, - Unknown 
Deflection: Үс = 0.0 Y, - Unknown 


By supplying these initial values to the program, we 
obtain the computer output shown on the next page. The last 
column shows the exact deflections. In these results, note 
that all the initial values are exact, and that the last value 
found, namely the deflection of the beam at x = 1.0, 1S expected 
to have the biggest cumulative error. However, as we can 


see, we get 3 to 4 digits accuracy. Here we can conclude 
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EX 2 "PE 


CLAMPUUSPRES ЖЕЙДІ 
BCUNDARY CONLITICNS: Y(0)= 0.0С00 
M(J)- 0.3159 

X: SFEAR Mev ery ЕР 
0.30 = Је о 20600 0.315000 0. COOQUU) 
0.05 —0.596000 6.350020 0.С18125 
9.10 = 225: 0008 C.325000 U.C35000 
0.15 -—0„506С600 С.300099 GC 20025 
9.2 С -0.596000 0.275000 0 ,C650230 
0.25 -06. 5CCCOQ С.250026 ТО 25 
0.30 —6- 565600 225520 0.C900200 
0.35 -0.5.6CQJ9 0.200000 0.100625 
9.40 о 0.175290 0.110490 
0.45 -0.520090 С.15У0929 О етен 
0.50 -2.525002 С.1252905 0.122909 
9.55 - 0.432 ССЈО 0.101250 0.130646 
0.00 -0.42C000 с .С 80009 0.135167 
0.65 Se ССОО 0.061259 0.138538 
0. 7С -0.3.CCO0 С.С42909 0.141334 
О. 75 -U 23060) 0.931250 0.143230 
О.С 46090 C.920900 0. 144501 
0.85 -9 130000 С.011250 0.145272 
0. 90 -0.1J3CCoQ9 С. ССБПОМ 0.145607 
О 20.020030 U. 201242 0.145313 
Е. ОС -0.02C000 -C.000001 0.145354 
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V(0)= 


(C)= 


LEFEL CEJ] 


C. 0000900 
C.J00433 
CULT 
C. 203937 
С 006556 
C. 01 0416 
С.Ј14625 
(91 9335 
C. 024660 
CU 90m» 
(. 036453 
С.04285 + 
С. 0+9 50% 
С. 250354 
6363951 
C.070475 
CeUTtoTOQ 
C.0384916 
2092 1I 
C.0994T9 
С.1267719 


De 2099 


-0.50С0 


EXACT DEF. 


-0.200000 


C.200458 
0.201 792 
0.303933 
0.396333 
0.910417 
0.314625 
0.019896 
0.024607 
Q4, 22:0 3105 
0.236459 
0.042355 
0.549505 
0.356 5 
0.963858 
0.079475 
0.077071 
0.234915 
0.092 198 
0.299479 
О. ТОСНО 


that previous values of all other columns are even more 
BN curate. 
3. Example 2:  Clamped-Roller Beam 

Бенет со Риасор Тисе а statically indeterminate 
beam loaded by a moment of 10 units of weight x unit of length, 
acting at the middle of the beam. As before we are going to 
determine the shear, moment, slope and deflection of the beam. 

ша Solution 

This example will emphasize three things: How to deal 
with a) concentrated forces, b) concentrated moments and c) the 
necessary initial conditions required for this particular 
Case. 

а Concentrated Forces 

We can conceive of a concentrated force as a dis- 

tributed force of high intensity distributed over a very small 
length of beam. Recall from the previous chapters that the 
slope of any approximating line at any point i is given by 
(2.9) where f(x) is the distributed load for the case of beams. 
In this equation the integral in the numerator represents the 
area under the f(x) curve. This area is actually the total 
load in the specified interval. By decreasing the interval 
while still holding the same inside area constant (same total 
load), the distributed force tends to a concentrated force of 
equal magnitude as the total load. In the limit it becomes a 
concentrated force acting on a mathematical point on the beam. 


For our purpose, in the integrating algorithm we simply let: 
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m. = 0 for х > X; 
m. = E FOE х = х 
СОТА ВЕ ИЕ fOr х = х 


where P is the concentrated force acting at point Xi. 
b. Concentrated Moments 
A concentrated moment can be treated in a similar 
way if we use the second order differential equation. However, 
in our case we need to decompose the given moment into a pair 
of concentrated forces of equal magnitude acting in opposite 
directions at equidistant points from the location of the 


given moment. From statics we know that: 
Moment = Force x Distance 


Here we can have several combinations of force 
times distance provided we keep this product a constant equal 
to the given concentrated moment. However as we will see 
later, accuracy is achieved only when we use very small dis- 
tances and consequently very large concentrated forces. In 
this way of moment decomposition, the usual behavior of a 
concentrated moment is achieved. For larger distances the 
couple is not an accurate representation of a concentrated 


moment and the results are not so accurate. 
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In the current example we use M - 10 and it is 


decomposed as follows: 


10 = 1000 x 0.01 


a pair of opposite forces of 1000 units acting 0.01 units of 
length apart from each other. In our problem we use 200 
intervals of 0.005 each, one concentrated force is located 

at x = 0.495, and the other opposite force of equal magnitude 
is placed at x = 0.505. So the loading function to be inte- 


grated is defined by 5 partial functions as follows: 


£1(x) = 0 > 
J£,(x)dx = -1000 x = 0.495 

fa (x) = 0 в 5 05 
ff,(x)dx = 1000 х = 0.505 

Е5(х) = 0 Вто 


с. Necessary Initial Conditions 
= со considers now the B.C,'s. For the present 


configuration we have: 
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Shear: E = Unknown M = Unknown 


Moment: Vou = Unknown ма a") 
О n 
Торе: ' = ! = 
SLO Y 0 Үл Unknown 
Deflection: Y = Ү О 
О n 


To start the algorithm we need to determine the initial 
shearing force Үс and the initial bending moment Yo- From 
eguations (3.13) to (3.16), after replacing the known values 


ме сев: 


ћ5 


e 
|| 


+ nny + Yu 
О О 


2 
Y! = боб, + Nohy" + nhy" 
n l l O O 
0 = hs +N hy" « N.h^v" 
O O O ili O 


This is a set of three simultaneous algebraic equations with 
three unknowns and we solve for the initial conditions. The 


results are: 


h(N OS. п5 1) 
о nN, - N Een 
l О 
S - N,S 
DES Ue 12 
ща ^. nN, -N (5.5) 
l О 
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All the terms involved in these expressions are 
wormed By пе relations (3-13) to (3.22), and they can be 
computed in advance since they are independent of boundary 
conditions. These initial values are only approximations as 
we know, but serve to match the given B.C.'s at the other 


extreme. The results are: 


ше 2 е Exact = 11.25 


Shear: ү"! 
О 


-].21219 Exact = -1.25 


Moment : үп 
О 


Running the program with these initial values, we 
get the results shown on the next page. The deflection appears 
accurate to about 2 percent for the more meaningful values of 
deflection. However, if we reduce the distance between the 
Opposite forces, we improve the accuracy. But if we increase 
it, the discrepancies grow enormously even for a very small 
increment, and it does not reflect a concentrated moment 
behavior. 

The slope shows the expected behavior. However, we 
should not expect high accuracy since we have started the 
algorithm with approximate starting points of shear and moment. 
At x = 1.0, we have Ya = -0.631. The exact is -0.625. Pre- 
ceding values of slope are expected to be better since here 
we start with an exact known initial point. Similarly for 


the moment we start with an approximation but get more accurate 


o3 


ов E 
СЕВЕРЕ Е 


BCUNDARY CONCITICNS: Y(0) 22.0 Y(1)=0.0 
О-о ао 

x: SHEAR МСМЕМТ SLOPE OZFLECTIJI EXACT ШЕШ 
0.00 Тоа епо G.C90000 C.000900 0.200000 
0.05 11.212190 -С.651584 -0.046594 -С.291282 -0.201327 
Оа li.212190 -С.С5С975  -0.C65158 -С.004190 -0.00+374 
0.15 ШЕ 20) C.4696034 -0.(55692 -С.097327 -0.907734 
0.20 11.212190 1.030243 . -0.C18195  -C.009290 | -0.009999 
0.25 EO „530853 0.C47333 -С.028677  -0.009766 
0.3C IE 2.151463 0.140890 . -C.004087 -0.905024 
0.35 D 2. плата. 0.262479 C. 005382 0.003828 
0.4C LI. l210 3.2126000 0.412998 C.022630 0.020000 
0.45 IS 3.833291 0.539747 C.047561 0.044297 
0.56 -$88.727500  -3.108097 0.164177 С. 081965 0.973125 
0.55 11.212190 -5.045486 . 504137 C.113698 0.110390 
0.6C 11.212190 -4.484873 0.265877 62132755 C.130000 
0.65 li.212190 с 4269 0.655649 (.140755 0.138359 
DIU 11.212190 -3.363659 -0.126550 C.136867 0.130875 
0.75 11.212190 -2.803050 -0.280717 Са (28670 0.120953 
0.35 11.212190 — -2.242440 — = 0.0400 855 С.111265 0.119000 
0.85 11.212190 1.681830 -0.504961 C.283354 0.087422 
0.90 11.212190  -1.121222 -0.415037 C. 061239 0.060625 
0.95 11.212190 -С.560013 -9.617084 C.2031319 0.931016 
1.00 11.212190 -90.000003 -0.631998 C. 000000 с.200с00 


54 


as we proceed, since we end up with the expected result of 
WO at xX = 120. 

In the case of the shear force, we see it is con- 
stant throughout the beam length as expected but it is shifted 
by an amount of 0.0378 from the exact. It is interesting to 
note that at the point of the loading we have a very strong 
shearing force. For a pure moment acting at this place, we 
shall neglect this result. However, if we do have an actual 
couple than a big shearing at this point should be expected. 


Note that the algebraic sum of shear forces gives: 


HEIDE 98937575). «= L000. 0 


which is indeed the couple acting at this point. 
5. Example 3: Pin-Roller Beam 
Refer to Fig. 5.lc. This is a statically determinate 
Simple supported beam subjected to a triangular type of load 
distributed over the central portion of the beam as indicated 
in the figure. We need to determine the shear, moment, slope 
and deflection of the beam. 
6. Solution 
First of all, some comments about this configuration 
are in order. This example was chosen in order to emphasize 
the ease with which we can switch from one kind of load dis- 
tribution to another. As explained at the beginning of this 
ШЕ, ІБ 15 вле Басе that the algorithm solves a. "new" differ- 


ential equation at every single step that allows us to deal 
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with different loads. Therefore it is possible to have a 

different kind of loading per subdivision. Another reason for 

selecting this problem is that here we deal with linear forcing 

functions and we are going to perform four successive integra- 

tions of a linear function. This will lead to a polynomial 

of fifth degree and a bigger cumulative error is to be expected. 
Proceding with the solution, the distributed load 


f(x) to be integrated is given by four different functions: 


f} (x) = 0 Qu bon 02 
fj(x) = 4х-1 хо 
£4(x) = -4х+3 О БОС ЕИО 
£, (x) = 0 О ЗАРА 


As indicated above, the computer program is such that 
it can switch from one kind of loading to the other at the 
specified point. Now, in order to start the algorithm we 
need to determine the initial points from the B.C.'s of a 


simple supported beam. These are: 


Shear: Y^ = Unknown ren = Unknown 
Moment : Y AEE S) Y" = 0 

О n 
Slope: Y = Unknown Үл = Unknown 
Deflection: Y = 0 M = (0) 

О n 
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Now, from the known relationships after we substitute in the 


given values, we get: 


О = hS, + ПАУ 


2 


0 = h^s + N h^y" + nny: 
О О О О 


This is a pair of simultaneous algebraic equations in two 


unknowns. We solve for Y^ and Үс: The results are: 
12222 
' = iid = 
Ше = (=) (№52 05) (5.6) 
5 
y" = 2 2 (5.7) 
O n 


Similarly, the computer supplies the values of the 


variables involved in these relations and these are as follows: 


Shear: Ye = -0.12499 Exact = 0.125 
Moment: Үс = 0.014972 Exact = 0.014974 
which shows a remarkable accuracy. Running the program with 


these initial values we obtain the results shown on the next 
page. A comparison between the approximate deflection and the 
exact reveals an error of about 2 percent at the point where 
the largest deflection takes place. In this example, we 


expected a larger error, but the algorithm appears to perform 
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EXAM ES 
PIN=RACLLER BEIN 


BOUNCARY CONCITICNS: Y(2)20.0 Ү141)-9.0 
` м(0)=0.0 4(1)20.C 

X: SHEAR MOMENT SLOPE DEFLECTION ET 
0.00 Бо падане C.C00000 C.C14972 C. 000000 0.000000 
0.05 391019459999 5: ы ONIS C.090746 0.2200; 
0.10 -0.124999 -6.012500 05014347 6.001476 0.001477 
0.15 7024990200050 ОСОБА C. 902175 0.002 Е 
0.2C 20.124999 00500 Ое o 6.002327 0.002328 
0.25 202124999 | о ] 2] 5759 ОЕ С. 003417 0.00 
0.30 a e E ае 0.C09348 C.003929 0.003920 
"PEL 2091942999 199 09998055 ONUS C.004347 0. 
0.40 -0.015999 -0.-041745 0695958 C.004653 — 0.200000 
0.45 ЕОР E INS 0.C02534 6.004349 0.002 625 
0.50 0.506000 —-0.052075 0.600000 6.004914 0.004 785 
0.55 0.C45CJ0 -С.050910 -0.C02584 С.004849 0.204706 
05916 053300000. 20400045 00901 6.504653 0.004 609 
0.05 0.125000 -с<.043030 —0.C07333 С.004347 0.004 322 
2 0.190000. —O2034515 9955900992. С.003929 0.003920 
0.75 0.12000 OO ое C.003417 0.003418 
0.3C 0125900. =0. ] КО C.092827 0.002328 
0.85 0.125000 0.017501 СІЗ = Googie 0.0021 5 
0.90 0.125000 —0.012500 ПЕТ ЛӘ 0.001477 
0.95 0125.00 02600250 M ОАА 52000126 0.0095 
1.00 0.125000 C.CQ0000 -0.С14972 —C.900000 0.000000 
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fairly well even in this case. Nevertheless, the error should 
be expected at the "free" (unknown) B.C.'s since the other 
extremes are fixed. Here we make use of the fact that the 
loading is symmetrical and the results should be symmetrical 
too. This is actually the case in the computer printout, so 
we can conclude that the results are indeed accurate. 
7. Example 4: Clamped-Clamped Beam 

Refer to Fig. 5.ld. This is a statically indeterminate 
beam clamped at both ends and loaded with a totally arbitrary 
discontinuous load known only as a table of values at discrete 
points as shown on the next page. We shall determine the shear, 
moment, Slope and deflection of the beam. 

Solution 

Before we proceed with this example, let us remind 
ourselves that we are not restricted to handle only exact- 
integrable forcing functions. We can deal with arbitrary 
functions as well, and this example intends to illustrate the 
procedure. In this problem, the computer program reads in 
data from a table of values instead of obtaining them by 
evaluating a given function. With this data the program per- 
forms trapezoidal integration in the usual way. The forcing 


function is defined by three other functions as: 


£1(x) = 0 0 < x < 0.3 
f, (x) = (from data deck) 45 х 50 7 
fa (x) = 0 Dc cc c Del 
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Unfortunately, no exact solution has been obtained for 
this case and we are not able to compare the computer results 
with the exact. Nevertheless some conclusions can be drawn. 


The B.C.'s for a clamped-clamped beam are as follows: 


Shear: No. = ВОК Пома Yt = Unknown 
Moment: Үс = Unknown | = Unknown 
Slope: Y' = 0 ү = 0 

O n 
Deflection: Y =. 0 Y _ 

О n 


To determine the initial conditions, we use the known 
relationships and substituting in the known values, we obtain: 
2 2 


ж. "i и 
О П 51 -- Nh Y + ППУ, 


= 3 3 сет 2 " 
O = h TT Non Шест Nh Үс 


This is a pair of simultaneous algebraic equations 


with two unknowns. Solving for the initial conditions we 


чес: 
п(М_5„ - nS.) 
У" = = nrn nn Села) 
E Ер О 
Se = 
ү" = 2 = (5.9) 
= 下 О 
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For this example the values supplied у спеншсотпршеен 


авес 


Shear: ме: = -0.948411 


Moment: = 0290507 


The computer printout is shown on the next page. Here 
we have another column to show the arbitrary loading. Note 
that we are using 100 intervals, but actually every fifth is 
prumbess 

Perhaps the only way to analyze these results would be 
to check if the B.C.'s have been met or not. The results show 
that they have. Another clue to assure some accuracy would be 
to look at the maximum deflection and the minimum slope. The 
loading was intentionally concentrated on the middle of the 
beam so we can expect the highest deflection and minimum slope 
in this neighborhood since we have similar B.C.'s at both 
ends. Ме сап say that this too checks. Consequently мен ке 


that the results are likely to be trusted. 
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EXAMPLE 4 


MESUPCSSCEAUPOC 3:4M 
И ВАБУ СОЛО Ties ҮШ Есер ¥(1)=0.9 
о Y'(1)=9.0 
х: HE X SE 
0.05 ыл (СО з сени 0.230501 
0.25 ШІ 0550411 оо 
0. 1C O HM 564 Uz 0155956 
0.15 co All 0.083245 
PC EE CL 0 .Са0825 
0.25 ЗІ 422 243555 
NOSE 22212 222502 
0.25 ЖІ 702270 0095001 
BSc ЕЗІ 2707 27137-57 
0.45 200 40:11. 0.065628 
9.56 ОО ОООО 522252011 ы ТО 
0. 25 оо (ОО mo 82 
C 6.320200 с. 
0. 65 acc ЕЕ Еа 
9.70 в OK 2.063138 
9.75 72515000 ЕСЕП 5615555 
OEC Um 150 Дали ВО отс 
0.45 2.610000 1.011580 22002 
TRUTS ЭБ оо 1.011530 0.139230 
3.95 SEO 1.611590 0.189629 
E. ONSE со 1.011530 o) SUIT 
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Up 


C. 200000 
С.212349 
БӘСІ уе 
(2023996 
Сеја о 
67627935 
оо» 
С. 022025 
C. 2169530 
6. 09494 
(С. 730614 


= 020 2082 
Е. Во? 
CUE LUT 
65120542 
z 1234252 
BED AO LO 
JE041 
“2013953 
о ко 


629990099 


ЕЕЕЕЕС ТІН 


0.900000 
©, 002638 
9.309994 
0.202059 
0.305344 
0.994132 
0.396 105 
022307339 
0.203330 
Os Io) 
0.99>224 
5292959 
Q . 208440 
0.207483 
По 0202 25 
0.09+ 870 
0.003454 
020021352 
ОЗ 
9901 219 
-C.000000 


C. THE COMPUTER PROGRAM 


FILES SLOPE wATFIV Al 


whe Ho he -Ағ -і. -Ат he he che ve بل ب‎ Ao یھ یھ‎ 
P- P P P 


4$ 4# 4+ 
4$ 4$ 4 
4$ 4$ 4 
4b tb 4 
45 484 
45 4 dt 
4+ 4$ 4$ 
3$ 4+ 4 
4$ 4$ 4$ 
4 4 4+ 
4$ 3$ dt 
4$ 4: 4$ 
4$ 4+ tt 
44 4 


~ 
= 
ate 
> 
~ 
> 


+ db + 


he Ае ае че ае А whe ~ ~ 
P- ~ ~~ 
же THE CCMPUTER FRCGR 
te e Ae awe whe whe whe whe یف‎ 
Aw m o а у а ле 


dt dh 4+ 
4b 4* 4+ 


~ 
~ 
~ 
~ 
~ 
~ 


THIS PROGRAM CALCULATES TRE SFEAR ING FCRCE. BENCING MOMENT SLOPE ANO 
CEFLECTIOM OF EEAMS BY PERFCFEMING FCUR SUCCESIVE INTEGRATICNS ON 
TRE LCADING FUNCTION. THE FCLLOWING CONFIGURATICNS ARE CCNKSICEREO: 

be “CANTILEVER (CLAMPEC-FFE=) BEAM 

e2e-CLAPP=0-ROLLER EE AM 

J -SIMPELE SUPCRHTEO (FIN-ROLLER) BEAN 

4. -CLAMPED-CLAMPED BEAN 
THE LOADING CAN BE ANY МЏУВЕЗ СЕ OCISTRUBUTEO ANC/CR CCNCENTRATEO 
FGRCES. CONCENTRATED МОМЕМТ< САМ ВЕ DECCMPCSED INTC A CCUPLE AT THE 
MINIMUN PGOGSSIELE CISTANCE ACTING AT THE PGINT CF APPLICATICN GF THE 
MCNENT o 


THE PROGRAM CCNSISTS ОҒ TEREE ROUTINES: 

1. SUBROUTINE SLOFE 

2. SUBROUTINE INCCM 

З. SUBROUTINE LCAC 

whICH ARS CEFINED IN THE FOLLCWING FARAGRAPRHS 


+ 
4 
+ 


жазы MAIN PROGRAM sese 


4% 


Be ho uA, o te Se «Ағ aha che Ao 
AM momo я 


CS IN NUMBER CF STEPS C 
LS SUSBROUT INE INCCM FOR 
L я 


S SUBROUTINE SLCFE ҒО 


ЕРТКЕО GY THE INTEGER VARIALBLE М 
INITIAL CCNOITICNS 
FI 


NAL RESULTS 


AOD 
> > П 
CC Pp 


4 
+ 
4 
tt 
4 
4 
4 


SUBROLTINE INCCM $$2222592322229222239225 $2 


+ 


и 


IN THE INTEGER VARIABLE K WHICH SPECIFIES THE KINO QF PROBLEM 
THE FOLLCwING COCE: 

CANTILZVER BHEAN (INITIAL CONOCITICNS MUST SE SPECIFIEC) 
CLAMPED ROLLER BEAM 

СКА МРЕО -СКАМРЕО ВЕАМ 

SIMPLE SLPPCRTEC &EAM ' 


€ D 
m (M 
нини о 


D UN 


I 
in 


DARAMENTEKS ARE ОЕРТМЕС AS FOLLOWS: 
CALL INCON(N.YC.YIC.Y20. Y3C) 

NUMBEH CF STEPS 

INITAL LEFLECTICN 

INITIAL SLOFE 

INITIAL ВЕМС:МС MOMENT 

INITIAL SHERING FCRCE 


2 «+ RARKAP 


YO 

vic 
үгс 
Y 3G 


THIS SUBROUTINE SUPPLIES THE INITIAL CGONCITICNS 8Y CALLIMNMG SLCPE。 
FCR A CANTILEVER EEAM THIS STATEMENT 15 МСТ ExECUTED 


аата а ата аа а а аа а аала аа ала аа ааа ааа а а а аа аа ааа а а ааа ааа. 


+ 


UD it 


SESS SUUQCLTINE SLOFE стт2тФФостФФтФФФ2:ФФФФФ2 
FORM TwO TASKS: l» CALCULATES ANO SUPPLY INCCM wITH THE NECESARY 
PARAMETERS ТС ОЕТЕКМІМЕ INITIAL CCNCITICNS 
2, EXECUTES THE SOLUTICN OF TRE PRCELEM 


m 4) 
L + 


THIS SUGROUTINE IS A THE CODING OF THE RELATIONSHIPS GIVEN BY 
EGCLAT IONS 3.12 ТО 3.12. ТРЕ SULMMATICN TER¥S ARE CONTAINEO ІК 
TRE COMMON STATEMENT 


THE CALLING PARAMETERS ARE CEFINECO AS FOLLOWS: 
CALL SLCPE(N.YO.YIC.YZ2. Y3C. J) 


N = SAME AS IN INCCM 
YC = ” 
VIG = 
vec = * 
үзсе” 
J = TAKE TWO VALUES; FINO INITIAL CONOITICNS ONLY 


= О 


2 $ TC 
J s TC EXECUTE THE SOLUTION 


wes SUBROUTINE LOAD 35533 


tt 
% 
+ 

| 


усу сусу УСУ СУ СУС С СУ СУ С Сз Су СУ СУ Сз СУ СУ СУ С 


PPLICS SLOPE wITH TRE NECES 
A 


SUPP FCRMATICN GBTAINEO FRCM THE 
LCACING FUNCTICNe ІТ САМ ЕС А N 


L FUNCTICN CR A DATA CECK。 
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Site. SLOPE wATFIV At 


C MUST GE WRITTEN SPECIFICALLY FOR EACH NEw PRCBLE*. WHEN DEALING 
С WITH CONCENTRATED LOACS.e THE VARIABLE RIAT IN RCUTINE SLOPE IS 
Е саст EQUAL TC THE CCNCENTRATED (| САО CIRECTLY 。 AT THE POINT CF 
ACTICN, 3 
C THe CALLING PARAMETERS ARE CEFINED AS FCUUCWS 
G CALL LOAD(E.H2) 
С 日 = CURRENT VALLE GF THe VARIASL= XK (LENGHT) GIVEN BY RCUTINE SLOPE 
С H2= CONTAINS THE VALUE CF TRE LOACING FUNCTION АТ x=8 
С 
се аана и лл ино оен л аа А 
C MAIN PROGRAM 
> 
REAL YO.YtiC.Y2C.Yv3Q 
INTEGER J, 
COMMON SO,.S1.4S2. SNO. SNI, RN. ^ 
ВЕАО( 5.10 JN 
10 ЕОВМАТ(1 3) 
Hzl./FLJAT(N-1) 
Әмсг( ЗАТ(Һ-!) 
СА! ІМСОМ(М»Ү0» Ү10.Ү20.Ү30) 
J=0 
CALL SLOPE(Ne VYOeYVYIO~Y 20 e¢Y20 -J ) 
STOP 
END 
C 
п 
SU3HOUTINE (САО(Е.Н2) . 
> 
REAL detHi2 
F (x у= ооо 9 
~2=F (6) 
RETURN 
ENO 
C 
и 
SUBROUTINE 1 ЧССм( метО,„У1С .” 20 „У3С) 
(а иа поели ан пе селе па ле ла а пене пе ле ле пр mE MA WE А 
ВЕА( УО,У:С,У20,У30 
INTZGER „К, 
COMMON 50,51, 52, <мМ0. < №1 „КА, н 
ЭЯЕАО( 5.10 ) К. ҮС, 10 о Ү2С, ҮЗС 
10 РОВМАТ( [1 4F10.4) 
ІЕ( (К. 206 С) Ок. (ко Ст, 4) ) СО ТС 60 
ТЕ (к.ЕО.1) СС ТС 50 
ШЕТ 
CALL SLOPE(N.YC. Y tO. Y2O. Y 204. J) 
1Е(к-3) 20, 30, 40 
CLO > SSS #2 
Е CLA MPED-RCLL ER 
Е Що 
20 ян! ТЕ (6.2<) 
25 FORMAT(*1*,.*EXAMPLE 2*//lX.*CLAMPCO-QCLLER EEAM* //LX,.*3QUNDARY CON 
SDITIONS: Ү(0)-0.0 Ү(1)-0.0%//27х.%Ү (0)=0.0 М(1):-0.0%//) 
ү2т-н2(5225М0-9М%650)/ (РЬКФ<М1-5М0) 
vage (S50-SN17S2)/( RN*SNMNI-SNO) 
RETURN 
ое 
c CLAMPSO-CLAMPED 
Е 
30 янт ТЕ (6.35) 
35 FORMAT(*1*,.,*—XAMPLZ 4*//1X. 'CLAMPZC-CLAMPED ВЕ АМ И/ ТХ, "ВОЦМОАКУ СО 
NOI TITONS: Y(0)=0.0 Ү(1)-0.0%//2<х,%Ү (0)-0.0 Ү (1):2020%//) 
Ү20-:нНн%(50%5М1-517%5М0)/(вм%5М0-5М15565М1) 
У ЗО = (St *SNE-R«4-7SOJ/Z(RNSSNO-SNITSNI) 
RETURN 
о 
ас же! ТЕ (6,4<) 
45 FURMAT(*1*,*€XAMPLF 3*//l1X4., *D IN-SQLLER BEAM* //1X. *SBCUNDARY CONOITI 
SONS: Ү(0)9.0 Ү(1)20.0%//27х.%М(0)-0.20 м(1)=0 0/4) 
ҮүІ0= ( (нием) х2) = ( 5мо= 52 - НМУ 50 ) 
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FILE? SLOPE wATFIV Al 


ҮЗО = —S2/fKRN 
RET URN 


Ae se uie е whe wate athe بی مض اا‎ «Ар «До Meats Ñe ate ante 
RE АД ee eee ee 


C 

C СКА МРЕО- ЕБЕЕ 
ле ЗА He ose ceste ue cho sp so ess ate te abe ale ate wha whe o 

Ст ap سیت‎ — Aue c Зе Тур ылы کہ سپ‎ 


SC ФНІТЕ(6.5<)ҮС.Ү1 

55 FORMAT(* 19. .*EXAM 

ZOITIONS: Y( 

-ах."У(О)%.Ғ8а.,4/ 
GO TO Zo 


4 
it 
it 
+ 
% 


x O DO +t 
m dh 
иш? 
№ 4 

o 


кл 
с 
2 
0 
о 
Е 
ч 
m 
z tt 
m + 


REAL SUJOs.SUOI.«SUO2. SUO3. SUII  SUI2 « 
*XRINTI«RINT2.RINT J4H1«-2. R8 [+ 65, SNMI SNM 
INFEGER Tede Keble NeMeNSIK eNSOL «NSIL«.IMI 
COMMON SQ»S1.« S2« SNO, SN1 eRACH 
5000292. 

500 1= 0, 

SUO 230. 

SUO 3-0. 

SUI 170. 

SUI 2-70. 

SUL 3=0. 

SU2 2z0. 

SU2 3-0. 

503309. 

RINT=O64 

ВІМТІ =0. 

RINT2zO. 

RINT 3-0. 

NS1K=0 | 
NSOL=0 

NSI1L 20 

H1=0.0 

H2-7-0.0 

IF( 4 GT.0O) GO TO 100 

ШЕТТЕ (6.2С) 

100 OC 400 i=leN 

1МЕ=1-1 

RI=FLOAT( IM1) 

взе {ISH 

CALL LOAD(B.H2) 

1iF(I1.«zQ.1)GC TO 200 
RINT=0.,S*+H“ (HL +H2) 

200 RINTA=RINTÍ 

RINTI=RINT2 

RINT2-zRINTI 

RINTI=R INT 


SU3J-ZSUJJ*PINTI 
S3 =SU33 


SU2 3=SU2 3+RINT2 
DOSE 5022 +5023 
22.5503 1+5022 


K=I-—2 

ІҒ( 1.16.2) К=0 

SUL 3=SUL3+RINT3 

501 2=501 2 +5013 
501125011 +5012 

NSI KZNSIK *K 

SNM L=FLOAT(NSIK) 

SNL -z2O0O.5*FI*SNMI 
S1-429575U33*SU22*SUI I1 
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ме //1X. * BCUNDARY CON 
Х4//27Х. *M(O)z*.FB8.^, 


eSU2J.SUJJ.RINT, 


FILE: SLOPE МАТЕ ІМ AL 


с 

шет-3 

IF( I-LE. 3) 120 

SUO JzSUOQ 3*RINTA 
5002-5002%5003 
500125001 »5002 
5000 = 5900 +5001 
NSIL-NSIL*L 

NSÜLZNSOL *NSIL 
SNM2-FLOJAT(NSOL) 

SNO z204.25*RÍI * SNMI*SNM2 
50-.125:5033%.75%55022%1.5%5011%5000 
IF( J« GT.,0O) со ТО 300 


Y32S3*Y iQ 
Y2-2H*T(S2*RISY30)*Y2O0 
YI=((0S1+SN1IF VIC) FHERIFY2C ) FH VIO 
Y-Z(((SO*SNO-*Y3Q0)TH*SNI* Y2C)7H*RITYAO) *H«* YC 
УНТТЕ (©. 1С) 6. УЗ. У. УТУ 
300 HIi=H2 
400 CONTINUE 
IF(J.GT.0) GO TC 999 
РІ ТЕ( 6,96) 
10 ҒОРВМАТ( 1Х.Ғ4.2.-2Х,аҒ12.6/) 
20 ҒОСМАТ( 2Х,%Х:%. ЭХ. SMEAR?” OX ,. *MCMENT*.6X,. *SLCEE* ,5X. *OEFLECTION' ) 
SG FOPRMAT(?*1*) 
99G RETURN 
ENO 
С 
SENTRY 


6/ 


VI. CONCLUSIONS AND RECOMMENDATIONS 


Based upon the research carried out in this thesis and 
the results obtained, the following conclusions can be drawn: 
1) First of all, the method, in the way it has been 

developed, shows the fact that the solutions are totally 
independent of each other. Any value of the unknown can be 
found without generating all previous ones. A closer look at 
equations (3.13) to (3.16) reveals that we can apply any of 
these relationships at any point i (O s нш атте т т 
vided we know the initial conditions and the summation terms 
which can be generated in advance. In all of these equations 
there is only one summation term which depends on the given 
function f(x), the other summations are series of integers 
totally independent of the given problem. This fact represents 
a good saving in terms of computer time if it is conveniently 
exploited. 

2) As it has been shown in the examples, the power of 
the method is perhaps its ability to deal with arbitrary func- 
tions and this is important since many engineering problems 
lead to these kinds of functions. 

3) So far, there has not been any error correction in a 
strict sense. The only corrective measure has been step size 
reduction. However, equation (4.9) shows that correction can 


be performed provided we know how the local error behaves. 
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When integrating constant and/or linear functions, (4.9) 
applies, but for other cases it does not. In any case, equation 
(4.3) indicates that the local error depends only on the second 
integral which in turn is one degree higher than the given 
function itself. This suggests the idea that if we know how 
the original function behaves, linearly, quadratically, etc., 
we could, to some extent, predict the error behavior and carry 
Out a correction. Further research is recommended in this 
particular case. 

Nevertheless, as shown in the examples, some accuracy has 
been achieved even working in single precision. For more 
complicated problems involving complex functions and requir- 
ing very accurate solutions, double precision is still a good 


possibility. 
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